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The statistical mechanics of phase transitions in dense systems of polydisperse particles presents distinctive 
challenges to computer simulation and analytical theory alike. The core difficulty, namely dealing correctly 
with particle size fractionation between coexisting phases, is set out in the context of a critique of previous 
simulation work on such systems. Specialized Monte Carlo simulation techniques and moment free energy 
method calculations, capable of treating fractionation exactly, are then described and deployed to study the 
fluid-solid transition of an assembly of repulsive spherical particles described by a top- hat "parent" distribution 
of particle sizes. The cloud curve delineating the solid-fluid coexistence region is mapped as a function of the 
degree of polydispersity 5, and the properties of the incipient "shadow" phases are presented. The coexistence 
region is found to shift to higher densities as 6 increases, but does not exhibit the sharp narrowing predicted 
by many theories and some simulations. 



I. INTRODUCTION 

A complex fluid is described as "polydisperse" when 
its constituent particles are not all identical, but exhibit 
a spread of size, shape or charge. Polydispersity is a per- 
vasive feature of both natural and synthetic macromolcc- 
ular systems such as colloids, polymers and liquid crys- 
tals. But for many years it was regarded as a practical 
and conceptual nuisance to be minimised in experiment 
and ignored by theory. More recently, however, it has 
become clear that polydispersity is actually a matter of 
fundamental interest and practical importance in its own 
right, giving rise as it does to a rich variety of phenom- 
ena not observed in monodisperse systems (for a review, 
see Sollich 1 ). Examples include novel phase behaviour 
such as an extreme sensitivity of phase boundaries to 
the presence of rare large particles, 2 critical points that 
occur below the maximum coexistence temperature, 3 ' 4 
density dependent wetting transitions 5 and surface size 
segregation effects, 6 ' 7 as well as interesting dynamical ef- 
fects such as an enhanced propensity to glass formation 8 
and the possibility of multistage relaxation processes. 9 

Despite the upsurge of interest in polydispersity in- 
duced phenomena in complex fluids, several basic ques- 
tions remain. A prime example is the nature of the freez- 
ing transition for simple spherical particles with a spread 
of diameters, which in quantitative terms is conveniently 
characterized by a parameter <5 measuring the standard 
deviation of the diameter distribution in units of its 
mean. Intuitively one expects that polydispersity should 
alter the location of the freezing curve and destabilize 
crystalline phases. 10 However, the character and extent of 
these alterations remain unclear despite extensive study 
by experiment, 11 ' 12 density functional theories, 13 ' 14 sim- 
plified analytical theories 10, 15-19 and simulation. 20-32 On 
the experimental side there is evidence for a "terminal" 
degree of polydispersity St above which a fluid will not 
crystallize. This is supported by a number of theoretical 
and simulation studies, some of which suggest that the 



terminal polydispersity arises from a progressive narrow- 
ing of the fluid-solid coexistence region with increasing 
5, with the boundaries of this region meeting at a point 
of equal concentration. 15,18 ' 26,31 ' 32 A reentrant melting 
scenario has also been proposed whereby compressing a 
crystal with a polydispersity slightly below 5t can cause 
it to melt. 18 One simulation study predicts the occur- 
rence of a partly crystalline "inhomogeneous phase" 27 at 
high polydispersity, while other work has argued that the 
suppression of crystallization is a dynamical effect arising 
variously from the low diffusivity of large particles, 33 the 
intervention of a glass transition, 30 ' 34,35 or anomalously 
large nucleation barriers. 36 

In our view, the reason for the diverse predictions of 
theory and simulation regarding the equilibrium phase 
behaviour is that much previous theoretical and simula- 
tion work has failed to cater properly for a key feature of 
polydisperse phase equilibria, namely fractionation - the 
phenomenon whereby the distribution of particle sizes 
can be different among coexisting phases, whilst the over- 
all distribution (across all phases) is fixed. One approach 
that does incorporate fractionation exactly (within the 
context of a mean field theory) is the moment free en- 
ergy (MFE) method. Previous MFE calculations by one 
of us 37 ' 38 find neither evidence for narrowing of the co- 
existence region with increasing 5, nor reentrant melting. 
Instead they predict that a polydisperse fluid can always 
split off a small amount of a solid phase having a nar- 
row distribution of particle sizes. The purpose of the 
present paper is to compare the findings of MFE calcu- 
lations with the results of tailored Monte Carlo (MC) 
simulations that similarly cater exactly for fractionation 
effects. 

The paper is arranged as follows. In Sec. II we outline 
principal aspects of the phenomenology of polydisperse 
phase equilibria. We then present in Sec. Ill a detailed 
discussion of general issues surrounding the best choice of 
simulation ensemble for obtaining accurate estimates of 
phase coexistence properties. Sec. IV describes the model 
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system we have studied and the bespoke techniques we 
have developed to determine fluid-solid coexistence prop- 
erties. The results of applying these techniques to poly- 
disperse spheres are described in Sec. V. Sec. VI sum- 
marizes our conclusions. 



II. PHENOMENOLOGY OF POLYDISPERSE PHASE 
EQUILIBRIA 

Typically one describes the polydispcrsity of a given 
system in terms of a density distribution //°' (a) which 
counts the number of particles per unit volume whose 
value of the polydisperse attribute a, lies in the range 
a . . . a + da. 39 In most real polydisperse systems the form 
of p^(a) is fixed by the synthesis of the fluid, and only 
its scale can change according to the degree of dilution 
of the system. Accordingly, one writes 

p^{a)=n^f{a) (1) 

where f(a) is a normalized fixed shape function and 

is the overall number density. Varying corresponds 

to scanning a "dilution line" of the system. 

At coexistence, particles with different a values will be 
partitioned unequally between the phases. This is the 
phenomenon of "fractionation". To describe it, it is nec- 
essary to define separate "daughter" density distributions 
p^\a) (7 = 1,2, . . .) which measure the distribution of 
the polydisperse attribute for each phase 7. When the 
polydispersity is fixed, particle conservation implies that 
the weighted average of the daughter distributions equals 
the fixed overall density distribution, or "parent" p(°\a), 
i.e. 

/)( a ) =n (»)/(a) = ^^V 7) W, (2) 
7 

with £( 7 ) the fractional volume of phase 7, where 
S 7 £^ = 1- This expression represents a generalisation 
of the lever rule to polydisperse systems. 

To illustrate the profound differences between phase 
behaviour in monodisperse and polydisperse systems, it 
is instructive to recall first the familiar case of the bin- 
odal curve of a monodisperse system in the density- 
temperature plane. This curve serves a dual purpose: 
on the one hand it describes the range of overall densi- 
ties for which phase coexistence occurs; and on the other 
hand it identifies the densities of the coexisting phases 
themselves. Now, for a polydisperse system the range of 
overall (parent) densities that leads to coexistence is sim- 
ilarly delineated by a curve in rj( '-temperature space - 
the so-called "cloud" curve. However, the densities of the 
coexisting phases themselves do not in general coincide 
with the cloud curve. Instead, as one varies the parent 
density through the coexistence region at a fixed 
temperature (say), one generates an infinite sequence of 
pairs of differently fractionated coexisting phases. 



Insight into this phenomenology can be gained by con- 
sidering the simplest case of a bidisperse (binary) mix- 
ture with densities p\ and pi of two species of particles 
of different sizes; these are the analog of p(°\a). Fig. 1 
sketches an isothermal cut through an exemplary bulk 
phase diagram, showing a region of fluid-fluid phase sep- 
aration with tie- lines that shrink to zero at a critical point 
(c.p.). 40 The dilution line constraint of a fixed shape for 
p(°)(<7) reduces in the bidisperse case to a fixed ratio 
Pi°Vp2^ (indicated by the dashed line in Fig. 1). As 
= plj -* + p^ is increased from zero, the system fol- 
lows the dilution line into the coexistence region, which 
it enters (and leaves) at a "cloud point" . For a given 
point on the dilution line inside the coexistence region, 
the parent splits into two daughter phases located at the 
ends of the tie-line which intersects this point. However, 
owing to fractionation, the daughters lie off the dilution 
line. 



P2 







shade 


w point 








1 / cloud point 








/ / // C ' P * 


/' cloud point \ 
/ shadow point 



Pi 

FIG. 1. Schematic of an isothermal cut through the fluid- 
fluid phase diagram of the binary mixture described in the 
text. Thick solid line: coexistence boundary; thin solid lines: 
tie-lines; dashed line: dilution line, here chosen as pj ' = p£ ■ 

Since an infinite number of points occupy the dilution 
line between the cloud points, an infinite sequence of dif- 
ferently fractionated coexisting states is encountered as 
the coexistence region is traversed. Customarily one sin- 
gles out the end points of this sequence for special at- 
tention, i.e. the case of incipient phase separation that 
occurs when the value of coincides with one of the 
two cloud points. Under these conditions, one daughter 
phase has a fractional volume of essentially unity and 
consequently (from the lever rule, Eq. (2)) a density dis- 
tribution that is identical to the parent, while the other 
phase - known as the "shadow" - has an infinitesimal 
fractional volume and a density distribution that deviates 
maximally from that of the parent. The curve formed by 
plotting the number density (or packing fraction) of the 
shadow phase as a function of temperature is known as 
the shadow curve. 

Qualitative differences between the phase behaviour 
of monodisperse and polydisperse systems are also evi- 
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dent in other projections of the full phase diagram. For 
instance, in the pressure-temperature plane, coexistence 
for a monodisperse system occurs along a line. By con- 
trast for a polydisperse system, coexistence occurs within 
a region of the pressure-temperature phase diagram. 41,42 
This is because each of the infinite sequence of coexisting 
states along the dilution lines is associated with a differ- 
ent pressure, and so the pressure varies as one crosses 
the coexistence region. For a fuller account of the phe- 
nomenology of phase behaviour in polydisperse fluids, we 
refer the reader to a previous review. 1 



III. SIMULATION STRATEGIES FOR POLYDISPERSE 
SYSTEMS 

Computer simulation is a principal route to deter- 
mining the phase behaviour of model thermal systems, 
and many such investigations of polydisperse systems 
have been reported. 20-32 However the results that have 
emerged from the various studies are generally much 
more at variance with one another (even for a given 
model) than is typically the case in comparable studies of 
monodisperse systems. It therefore seems appropriate to 
consider carefully what particular issues and challenges 
polydispersity presents for simulation and how best to 
tackle them in practice. To this end we provide here 
a detailed critical assessment of the relative utility of 
the various statistical ensembles and simulation strate- 
gies that have been employed previously for the compu- 
tational study of phase transitions in polydisperse sys- 
tems. Our discussion distinguishes two categories of en- 
semble: those that fully permit density fluctuations at 
the level of individual particle species ("unconstrained 
density ensembles" ) , and those ( "constrained density en- 
sembles") that do not. We argue that the latter category 
are seriously deficient both with regard to their ability 
to probe equilibrium phase behaviour, and their suscep- 
tibility to finite size effects. By contrast, unconstrained 
density ensembles - whilst requiring some effort to fully 
realize their benefits - constitute the method of choice 
for accurate studies of polydisperse phase equilibria. 



A. Unconstrained densities 

As is now well established, 43 ' 44 an efficient and ac- 
curate strategy for determining phase boundaries in 
monodisperse systems is to employ an ensemble that cap- 
tures on a global scale the fluctuations characteristic of 
the transition. Prime examples are the grand canonical 
ensemble (constant fj,, V, T) and the isobaric-isothermal 
ensemble (constant N,p,T). 45 The common feature of 
these ensembles is that they avoid physical contact (i.e. 
interfaces) between coexisting phases even when operat- 
ing at a coexistence state point. Instead the two phases 
are linked - and their relative free energies measured - 
via a phase space path 46 (the negotiation of which may 



entail biased sampling) that permits the simulation to 
fluctuate back and forth between the disjoint configu- 
ration spaces of the pure phase states. Although this 
path may traverse mixed phase (i.e. interfacial) states 
en-route between the pure phases, the associated surface 
tension penalty is such that at any one time, the sys- 
tem will be found with overwhelming probability in pure 
phase states. Accordingly there is ample opportunity to 
sample the properties of the phases whilst they occupy 
the maximum available system volume - thus minimiz- 
ing finite-size effects. Indeed it has been shown that when 
the system is found with equal probability in each phase 
(the "equal weight criterion"), finite-size corrections to 
measured coexistence properties are exponentially small 
in the system size. 47 

In seeking to study phase behaviour in polydisperse 
systems, it is clearly desirable to retain the aforemen- 
tioned benefits of the (p,V,T) and (N,p,T) ensembles, 
whilst at the same time generalizing them to permit sam- 
pling of the fluctuations in the polydisperse attribute. 
More specifically, one wishes to cater for fractionation so 
that the distribution of a (which for definiteness we shall 
take as the particle size) can vary from phase to phase. 
For the particular case of MC simulations within the 
(/U, V, T) ensemble, the total particle number fluctuates 
(by means of insertions and deletions) and incorporating 
polydispersity entails imposing a distribution of chemical 
potentials /x(tr) and additionally introducing particle re- 
sizing updates. Doing so permits fluctuations not only in 
the overall instantaneous size distribution p(o~), but also 
in the daughter distributions, thus facilitating efficient 
relaxations to a fractionated state. 

However, to tackle the experimentally relevant scenario 
of fixed polydispersity, one further needs to ensure that 
the ensemble averaged density distribution (across all co- 
existing phases) corresponds to some prescribed parent 
as expressed in Eq. (2). To do so entails solving an in- 
verse problem for p(a), 48,49 a task which - at first sight - 
appears complicated by the fact that each phase will (in 
general) either be absent or occupy the whole simulation 
box rather than its true canonical fractional volume 
Fortunately, a straightforward solution has been devel- 
oped by Buzzacchi et al. in which /it(cr) and the ^ 7 - ) are 
treated as parameters which are tuned iteratively such 
as to simultaneously satisfy an equal weight criterion for 
the coexisting phases and the lever rule constraint. 50 This 
technique permits the determination of coexistence prop- 
erties with finite size effects that are exponentially small 
in the box size, even in the limit in which one of the 
phases has an infinitesimal fractional volume. 

Whilst appropriate for studying fluid phase transitions 
in polydisperse systems, the (^(c), V,T) ensemble is un- 
suitable in situations where one or more of the coexisting 
phases is a crystalline solid because the number density 
of a crystal cannot easily fluctuate without engendering 
defects. Instead it is advantageous to appeal to a hybrid 
of the grand canonical and (N, p, T) ensembles known as 
the isobaric-semi-grand canonical ensemble (SGCE). 45,51 
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This is the analog of a monodisperse (N, p, T) ensemble 
where the particle number is fixed but the prevalence of 
the various particle sizes is controlled by imposing chem- 
ical potential differences ft(a) that are measured relative 
to the chemical potential of some reference particle size. 
The SGCE has been widely deployed in the simulation of 
fluid mixtures 45 and was first applied to polydisperse sys- 
tems by Kofke and Glandt. 51 Although the overall par- 
ticle number is fixed in the SGCE, via the constraint 
V J p(a)da = N, particle size updates nevertheless per- 
mit the sampling of many realizations of the polydisperse 
disorder, while updates of the overall volume allow the 
total system number density to relax in each phase. Thus 
the ensemble provides as many degrees of freedom vis a 
vis density fluctuations as the {n(cr), V, T) ensemble and 
hence permits efficient treatment of fractionation. An ad- 
ditional attractive feature of the SGCE is that it provides 
direct access to the coexistence pressure which is not the 
case in the (fi(a),V,T) ensemble. In Sec. IV D we shall 
describe how to combine the SGCE with the method of 
Buzzacchi et al. 50 to permit accurate determination of 
coexistence properties of polydisperse systems. 



B. Constrained densities 

We now consider the comparative utility of other en- 
sembles that have been utilized in the study of poly- 
disperse phase equilibria, but which constrain the par- 
ticle densities to a greater or lesser degree. These are 
the microcanonical (N, V, E) 30 canonical (N, V, T) 27 and 
isobaric-isothermal (N,p,T) 31 ensembles. We focus first 
on the fully constrained (N, V, T) and (N, V, E) ensem- 
bles. Here two or more extensive quantities are con- 
served and consequently at coexistence the simulation 
box divides into separate regions, each occupied by one 
phase which is connected by an interface to its neigh- 
bour. As a result, measurements of the properties of one 
phase can be affected by the presence of the interface 
with the other phase. Whilst in a monodisperse system 
- where the properties of the two coexisting phases are 
independent of their fractional volumes in the thermo- 
dynamic limit - one can mitigate interfacial effects by 
measuring the properties at equal fractional volumes, i.e. 
£(i) _ £(2) _ 2/2, the situation is much more delicate for 
polydispersity. Here, owing to fractionation, the charac- 
ter of the phases depends inherently on their fractional 
volume, and accordingly their properties need to be de- 
termined even when one of the phases occupies a small 
fractional volume. However, in a finite-sized system, this 
generally translates to a small absolute volume, and thus 
the properties of such a phase will unavoidably be dom- 
inated by its interface. This effectively precludes the ac- 
curate determination of bulk properties in these ensem- 
bles. In particular, near cloud points the corresponding 
shadow phases may occupy too small a volume to form 
at all, because of interfacial free energy costs, and this 
can lead to serious misestimates of cloud point locations. 



Finite size effects are further exacerbated in the 
(N,V,E), (N,V,T) and indeed the (N,p,T) ensembles 
by the fact that polydispersity is incorporated by assign- 
ing each of the N particles a fixed size drawn from the 
prescribed parent distribution. Consequently only a sin- 
gle realization of the parent is considered rather than a 
fluctuating sample as occurs in the {p(cr), V, T) or SGCE 
approaches. Fixed particle sizes represent a particular 
handicap in the interface-forming (N, V, E) and (N, V, T) 
ensembles when fractionation effects are strong, so that 
e.g. one daughter phase has a distribution that is strongly 
peaked towards the largest particles. 38 ' 52 If in that parti- 
cle size range the parent density is low, then the lever rule 
forces the relevant daughter phase to have a small frac- 
tional volume and hence also small absolute volume, with 
the problematic consequences discussed above. Equally 
if not more importantly, the sampling of the size distri- 
bution with N particles of fixed size may be too coarse 
to represent such sharply peaked daughter distributions 
accurately, especially in a tail region of the parent. This 
causes deviations from thermodynamic bulk behaviour 
that cannot be corrected in any straightforward way. 

A further practical disadvantage of fixed particle sizes 
is that for dense systems where diffusion is inhibited, par- 
ticles of a given size may not reach their favored phase 
on simulation timescales. Whilst in MC implementations 
this can be overcome by employing long ranged particle 
moves or particle exchanges, 27 it effectively prevents re- 
laxation in molecular dynamics (MD) simulations which 
employ realistic dynamics. A case in point is the re- 
cent work of Nogawa et al., 31 who performed MD in an 
(N,p,T) ensemble to simulate a system of size-disperse 
spheres initialized in an interfacial state of coexisting 
fluid and solid phases. As discussed above, interfacial 
configurations are disfavored over pure phase states in 
the (N, p, T) ensemble. But once initialized, an interface 
can nevertheless be maintained for a substantial period 
if the prescribed pressure is carefully tuned such that the 
interface is stationary (neither grows nor shrinks on aver- 
age) . This pressure serves - in principle - as a measure of 
the coexistence pressure. However, because in the study 
of Nogawa et al. the size distribution was initialized to be 
identical in both phases, and since no significant fraction- 
ation of the phases occurred on simulation timescales, the 
results emerging from this study cannot be regarded as 
representative of equilibrium. 

Further evidence for the disadvantages of employing 
constrained density ensembles is to be found in the re- 
cent work of Fernandez and co-workers 27,30 who em- 
ployed (A, V, E) and (N, V, T) MC simulations to study 
the fluid-solid coexistence of polydisperse soft spheres. In 
these ensembles, an interface between coexisting phases 
is mandated to form in the thermodynamic limit, as 
pointed out above. Apparently, however, none was ob- 
served at small degrees of polydispersity, with the sys- 
tem simply fluctuating between the pure phases - a fea- 
ture which presumably reflects the rather small system 
size used. Only at high degrees of polydispersity did an 
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interface form, but the authors interpreted its appear- 
ance as a polydispersity induced "inhomogencous phase" 
rather than an essential feature of their choice of ensem- 
ble. Even in this region of apparent phase separation, 
no account was taken of the finite width of the coexis- 
tence region, i.e. no attempt was made to distinguish the 
infinity of coexistence curves, or even the boundary of 
the coexistence region, either in density or temperature 
space. This presumably reflects the difficulties of deal- 
ing correctly with situations where one phase is incipient, 
as described above. Although use of long range particle 
exchange moves resulted in some evidence for fraction- 
ation well inside the coexistence region, it is our view 
that the results emerging from Refs. 27,30 are nevertheless 
qualitatively unreliable. Indeed it was this conviction 
that catalysed the present study, in which we revisit the 
phase diagram of the same model studied by Fernandez 
and co-workers, but using the SGCE and specialized data 
analysis techniques in order to extract the correct phase 
behaviour. 



C. Fixed versus variable parent 

As discussed in Sec. II, in many complex fluids such 
as colloids and polymers, the form of the polydisper- 
sity is determined by the process of chemical synthesis 
and hence the shape of the parent density distribution is 
fixed. The phase behaviour is then as described in Sec. II. 
However, most previous computational studies that use 
the SGCE or (/x, V, T) ensemble to facilitate fractiona- 
tion have not sought to adapt the form of /i(cr) in order 
to ensure that the ensemble averaged density distribution 
p(a) had a fixed functional form, corresponding to a pre- 
scribed parental size distribution /(o-). 22 ~ 24 > 28 > 29 Instead 
the activity distribution exp[/3//(cr)] was assigned a fixed 
form such as a Gaussian, peaked at some Co, and various 
widths of the Gaussian were used in order to change the 
degree of polydispersity. 

In such a fixed chemical potential approach, p{a) can 
vary dramatically across the phase diagram with the re- 
sult that one doesn't capture the phase behaviour of a 
particular parent form as would be the case experimen- 
tally. Moreover, many of the characteristic features of 
phase coexistence for a fixed parent are absent. Specifi- 
cally, when crossing a coexistence region at fixed chemical 
potentials, the system follows a tie line rather than cut- 
ting an infinity of tie lines as occurs when the dilution 
line constraint (Eq. 1) is imposed. This difference is then 
manifest in the fact that coexistence occurs only along 
a line in the (T,p) plane for fixed chemical potentials, 
rather than within a region. Such features are more akin 
to those occuring in a monodisperse system and accord- 
ingly the fixed chemical potential approach misses much 
of the essential phenomenology of experimental systems, 
a fact that severely limits its applicability. 



D. Summary 

We conclude this section by distilling the salient points 
of the above commentary. The central message is that in 
computational studies of phase behaviour in dense poly- 
disperse systems, the choice of simulation ensemble can 
have far greater implications for the severity of finite-size 
effects and the pace of relaxation than in monodisperse 
systems. Specifically, we believe that despite the fact that 
they ostensibly offer a simple route to fixing the parent 
form, polydisperse versions of the (N, V, E) and (N, V, T) 
ensembles are to be avoided. This is because on the one 
hand they are intrinsically hamstrung by serious finite- 
size effects (whose origin lies in interface formation) , and 
on the other hand they can suffer from very long relax- 
ation times. 

Suitable ensembles for dealing with a fixed parent are 
the (/i(cr), V, T) ensemble and the SGCE, the latter being 
most appropriate when solid phases are involved. Their 
strengths are two-fold, namely that they: (i) permit es- 
timates of coexistence properties with finite-size correc- 
tions that are exponentially small in the system size, even 
in the limit where one of the phases has an infinitesimal 
fractional volume; and (ii) accelerate fractionation via 
particle size updates that circumvent the physical diffu- 
sion process. The (modest) price to be paid in order to 
realize these benefits is the need to determine at each 
state point of interest a chemical potential distribution 
and values for the fractional phase volumes These 
must simultaneously satisfy the lever rule and an equal 
weight criterion. In Sec. IV D we describe how this can 
be efficiently achieved within the context of the SGCE. 

IV. MODELS AND METHODOLOGIES 
A. Models 

The systems that we shall consider in this work are as- 
semblies of spheres interacting either by a repulsive soft 
sphere potential (as considered by simulation) or a hard 
sphere potential (as studied in our moment free energy 
method calculations). The soft sphere interaction poten- 
tial between two particles i and j with position vectors 
Ti and Tj and diameters Oi and aj is given by 

v{r ij )=e{(j ij /rij) 12 , (3) 

with particle separation = |rj — r\,-| and interaction ra- 
dius <Jij = (<Ji + <Tj)/2. The choice of this potential rather 
than hard spheres is made on pragmatic grounds; in our 
isobaric SGCE simulations (to be reported below), any 
MC contraction of the simulation box that leads to an 
infinitesimal overlap of two hard spheres will always be 
rejected, so (particularly at high densities) we can expect 
higher MC acceptance rates using this "softer" potential. 
In common with hard spheres, the monodisperse version 
of our model freezes into an fee crystalline structure, 53-55 
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and temperature only plays the role of a scale: the ther- 
modynamic state depends not on and T separately 
but only on the combination n^^e/k^T) 1 ^. Phase dia- 
grams for different T then scale exactly onto one another, 
and we can fix e/ksT = 1. 

In all cases we consider parent size distributions of the 
top-hat form: 

(a)= ((2c)- 1 ifl-c<a<l + c (4) 
■' ^ ' 1 otherwise ^ ' 

Here the width parameter c controls the degree of poly- 
dispersity S = c/y/3, and we have set the mean particle 
diameter to 1. With these choices, and the interaction 
potential (3), our results are directly comparable to the 
phase diagram of Fernandez et al., 27 bearing in mind that 
in the latter work neither fractionation nor, at a more 
basic level, the presence of coexistence regions of finite 
width was allowed for. 



B. The isobaric semi-grand canonical ensemble 

Our simulations operate within the SGCE, wherein the 
particle number N, pressure p, temperature T, and a 
distribution of chemical potential differences fj,(a) are all 
prescribed, while the system volume V, the energy, and 
the form of the instantaneous density distribution p(o~) 
all fluctuate. 56 As discussed above, this is important to 
allow for separation into differently fractionated phases. 

Operationally, the sole difference between the iso- 
baric, semi-grand-canonical ensemble and the (JV, p, T) 
ensemble 45 is that one implements MC updates that se- 
lect a particle at random and attempt to change its diam- 
eter a by a random amount drawn from a zero-mean uni- 
form distribution. This proposal is accepted or rejected 
with a Metropolis probability controlled by the change 
in the internal energy and the chemical potential: 51 

p acc = min [1, cxp (-/3[A$ + {1(a) - , 

where A<E> is the internal energy change associated with 
the resizing operation. 

C. Phase Switch Monte Carlo 

Phase Switch Monte Carlo (PSMC) is a general 
method for determining phase boundaries but is particu- 
larly useful for dealing with fluid-solid coexistence. 57 The 
basic idea is to employ a reversible phase space leap that 
connects the configuration space of the fluid to that of 
the solid. This allows sampling of the disjoint configura- 
tion spaces of the two phases in a single simulation run 
and hence direct estimation of their relative free ener- 
gies. The method has been previously described in detail 
in the context of monodisperse systems, 58 as has its ex- 
tension to polydisperse systems. 59 We therefore confine 
ourselves to providing only a bare outline here. 



The method is based on a mapping between two ref- 
erence configurations - one for the fluid and one for the 
solid. A reference configuration for a given phase 7 is 
simply an arbitrarily chosen configuration of that phase 
defined by the associated set of particle sites {R} 1 ■ One 
can simply express the coordinates of each particle in 
phase 7 in terms of the displacement from its reference 
site, i.e. 

r] = R] + Ui . (5) 

Now, for displacement vectors that are sufficiently small 
in magnitude, one can clearly reversibly map any config- 
uration {r 7 } of phase 7 onto a configuration of another 
phase 7' simply by switching the set of reference sites 
{.R} 7 — > {i?} 7 , while holding the set of displacements 
{u} constant. This switch, which forms the heart of the 
method, can be incorporated into a global MC move. 

A complication arises, however, because the displace- 
ments {it} typical for phase 7 will not, in general, be 
typical for phase 7'. Thus the switch operation will 
mainly propose high energy configurations of phase 7' 
which arc unlikely to be accepted as a Metropolis up- 
date. This problem can be circumvented by employing 
biased sampling techniques to seek out those displace- 
ments {it} for which the switch operation is energeti- 
cally favorable. These are the so-called gateway con- 
figurations, which typically correspond to displacement 
vectors that are small in magnitude. 

Using the above formalism one can construct a sam- 
pling scheme which explores the configurations of high 
statistical weight in each phase, whilst returning at reg- 
ular intervals to the low weight gateway configurations 
that permit a switch to the other phase. In this manner 
one can directly compare the relative statistical weights 
of the two phases and hence determine free energy differ- 
ences and coexistence points. 



D. Fixing the parent distribution across coexisting phases 
and determining fractional volumes 

For SGCE simulations of a polydisperse system at 
some given N and T, we seek the pressure p and dis- 
tribution of chemical potential differences jj,(a) such that 
a suitably defined ensemble-averaged density distribution 
matches the prescribed parent p(°\a) = n^f(a). Unfor- 
tunately, the task of determining the requisite p and p,(a) 
is complicated by the fact they are unknown Junctionals 
of the parent. 49 To deal with this problem one can em- 
ploy a version of a scheme originally proposed in the con- 
text of grand canonical ensemble studies of polydisperse 
phase coexistence 50 and later extended to the SGCE, 59 
the latter implementation of which we now summarize. 

The strategy is as follows. For a given choice of 
and temperature T, one tunes p, p,(a) and the £^ itera- 
tively within a histogram reweighting (HR) framework, 60 
such as to simultaneously satisfy both a generalized lever 
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rule and equality of the probabilities of occurrence of the 
phases, i.e. 



/ (0 V) = E^ (7 V 7 V), 



£ = 0. 



(6a) 
(6b) 



In the first of these constraints, Eq. (6a), the ensem- 
ble averaged daughter density distributions p^(a) are 
assigned by averaging only over configurations belonging 
to the respective phase. The deviation of the weighted 
sum of the daughter distributions /5(a) = X) 7 ^ 7 V 7 H CT ) 
from the target n^°\f (°\cr) is conveniently quantified by 
a 'cost' value: 



A = 



p(a) - n (0) / ( °V) \ da . 



In the second constraint, Eq. (6b), 



,(7) 



n 



(7) 



(8) 



provides a measure of the extent to which the probabil- 
ity of each phase occuring, p^\ is equal for each of the 
n coexisting phases. Imposing this equality ensures that 
finite-size errors in coexistence parameters are exponen- 
tially small in the system volume. 47,50 

The iterative determination of p, p(a) and such as 
to satisfy Eqs. (6a) and (6b) proceeds thus: 

1. Guess initial values of the fractional volumes 
corresponding to the chosen value of . Usually if 
one starts near a cloud point, the fractional volume 
of the incipient phase will be close to zero. 

2. Tune the pressure p (within the HR scheme) such 
as to minimize A. 

3. Similarly tune p(a) (within the HR scheme) such 
as to minimize A. 

4. Measure the corresponding value of £. 

5. if £ < tolerance, finish, otherwise vary (within 
the HR scheme) and repeat from step 2. 

The minimization of A with respect to variations in 
p (step 2) can easily be automated using standard 1- 
dimensional minimization algorithms such as the "Brent" 
routine described in Numerical Recipes. 61 Similarly we 
used the "Powell" routine for the multi-parameter min- 
imization of £ with respect to variations in in step 
5. 62 In step 3 the minimization of A with respect to vari- 
ations in jl(a) is most readily achieved 48 using the fol- 
lowing simple iterative scheme for jl(a): 



ftk+i(c) = Pk{o) +a\n 



AO) 



(9) 



for iteration k — > k + 1. This update is applied simulta- 
neously to all entries in the histogram of fi(a), and there- 
after the distribution is shifted so that p(ao) = 0, where 
a o is the chosen reference size. The quantity < a < 1 
appearing in Eq. (9) is a damping factor, the value of 
which may be tuned to optimize the rate of convergence. 

The values of £^ and p resulting from the applica- 
tion of the above procedure are the desired fractional 
volumes and pressure corresponding to the nominated 
value of n(°). As has been described previously, 59 daugh- 
ter phase properties are obtainable by monitoring, sep- 
arately for each phase, the density distribution p^ (a) 
and the distribution of the fluctuations in the overall 
number density, P(n), and the volume fraction, P(rj). 
Here the volume fraction is defined in the obvious way as 
r] = J dap(a)(n /6)ct 3 , for a phase with density distribu- 
tion p(a). 



E. Analytical calculations: the moment free energy 
method 



Calculating analytically the phase behaviour of poly- 
disperse systems is a challenging problem. 1 This is be- 
cause for each of the infinitely many different particle 
sizes a one has a separate conserved density p(a). Effec- 
tively one thus has to study the thermodynamics of an 
infinite mixture, where e.g. from the Gibbs rule there is 
no upper limit on the number of phases that can occur. 

The moment free energy (MFE) method 1 ' 63 ~ 65 is de- 
signed to get around this issue by effectively projecting 
the infinite mixture problem down to that of a finite mix- 
ture of "quasi-species" . This is possible when the free 
energy density has a so-called truncatable form, 

/ = k B T J dap(a) {\n(p(a)) - 1] + / ex ({p0) , (10) 

where the excess part / cx depends on a finite number of 
moments of the density distribution, 

Pi = J dap(a)w l (a) . (11) 

This truncatable structure obtains for a large number of 
models of mean field type. Importantly for our purposes, 
it is also found in accurate free energy expressions for 
polydisperse hard spheres, with the simple weight func- 
tions u>i(a) = a 1 (i = 0,1,2,3). Specifically, we use for 
the fluid the Boublik, Mansoori, Carnahan, Starling and 
Leland (BMCSL) expression 39,66 ' 67 and for solid phases 
the free energy developed by Bartlctt 68 on the basis of 
the simulation data of Kranendonk et al. 69 for binary 
mixtures. 

The MFE method is a way of expressing the ideal con- 
tribution to the free energy from Eq. (10), which de- 
pends on the complete shape of the density distribution, 
in terms of the moment densities pi. The result is the 
moment free energy. The key feature of the method is 
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that if one then treats the quasi-species densities pi as if 
they were densities of ordinary particle species, and cal- 
culates phase equilibria accordingly, the results for cloud 
and shadow points are fully exact. The MFE approach 
is therefore the method of choice for our current investi- 
gation. We do not give further details of the numerical 
implementation here as these are set out in full in Rcf. 38 



V. FLUID-SOLID COEXISTENCE 

We have applied the simulation methodologies outlined 
in Sees. IVB-IVD to determine the fluid-solid coexis- 
tence properties of the polydisperse soft sphere model. 
In parallel we have employed the MFE method outlined 
in Sec. IV E to determine the coexistence properties of 
polydisperse hard spheres. The reason for not using soft 
spheres also in the analytical calculations is that there 
are no sufficiently accurate free energy expressions avail- 
able for this case. Qualitatively, however, we expect the 
hard and soft sphere systems to show similar behaviour, 
and will see that this is indeed the case. 

Our simulations of the soft sphere model comprised 
systems of N — 256 particles. The interparticle potential 
was truncated at half the box size and periodic bound- 
ary conditions were applied. Since for the system sizes 
studied the value of the potential is extremely small at 
the typical cutoff radius, no correction was applied for 
the truncation. As a preliminary step we determined 
the coexistence parameters of the transition from fluid 
(F) to face-centred-cubic (fee) solid (S) in the monodis- 
perse limit. This was done using the standard formula- 
tion of PSMC 58 with the results: 55 p cocx = 22.32(3), p F = 
1.148(9), ps = 1.190(9). Thereafter we attempted to lo- 
cate the fluid phase cloud point (£ = 0) for a narrow top- 
hat parent having c = 0.01. To this end we initialized the 
chemical potential difference distribution as /2(cr) = (for 
0.99 < a < 1.01) and p(a) = —100 otherwise, and as- 
signed the pressure the value p = 22.32 pertaining to the 
monodisperse limit. We then performed a long PSMC 
run, the results of which were reweighted (using the pro- 
cedure described in Sec. IV D together with the monodis- 
perse value n(°) = 1.148 as the initial guess for the fluid 
cloud point density), to yield accurate estimates of the 
fluid phase cloud point pressure, parent density and 
chemical potential difference distribution p(a), as well as 
the shadow phase daughter distribution. 

In order to progress to higher degrees of polydispersity, 
we proceeded in a stepwise fashion. The form of jj(a) pre- 
viously determined for c = 0.01 was extrapolated via a 
quadratic fit to cover the range 0.98 < a < 1.02, while 
the pressure for c = 0.02 was estimated by linearly ex- 
trapolating the results for c = 0.00 and c = 0.01. A 
new PSMC run was then performed, the results of which 
were reweighted to give accurate estimates of the cloud 
point parameters for the top hat parent with c = 0.02. 
In this manner we were able to steadily increase the de- 
gree of polydispersity, measuring the fluid cloud point 



pressure and density as well as the corresponding solid 
shadow properties. In an identical fashion we determined 
the dependence on polydispersity of the solid cloud point 
parameters, by finding the value of n^ ' for which £ = 1. 

For values of 5 > 0.069 on the solid cloud curve and 
8 > 0.087 on the fluid cloud curve, the system was found 
to make spontaneous transitions between the coexisting 
phases. Such an occurrence suggests a polydispersity in- 
duced lowering of the fluid-solid surface tension which 
normally provides the free energy barrier maintaining a 
simulation in the phase in which it is initiated. Unfortu- 
nately such spontaneous transitions undermine the phase 
switch strategy whereby one keeps track of the phase 
via the switch operation. Thus we were unable to reach 
larger values of the polydispersity than those quoted, but 
note that the values we can access are well within the 
range of typical experimental systems. 

In seeking to compare the simulation and MFE results 
for phase coexistence properties, it is necessary to adopt 
a means of translating values of S between soft and hard 
spheres. As standard mappings from soft to hard spheres 
fail at our densities, we use phase diagram topology to 
identify comparable points. This is achieved by scaling 5 
in order to match the location of certain phase transitions 
that occur in the solid region of the phase diagrams of 
our models 70 as will be described in detail elsewhere. 71 
The scaling is linear and maps S = 0.088 for soft spheres 
to 5 = 0.0546 for hard spheres. 

Our results for the cloud and shadow curves in the vol- 
ume fraction-polydispersity plane are shown in Fig. 2. In 
this representation one sees a close accord between cloud 
and shadow curves. This is surprising given the presence 
of fractionation. But one can confirm by perturbation 
theory 72,73 that the behaviour we observe comes from 
the fact that in our models polydispersity enters only via 
the particle size and not the strength of the interaction. 
For the MFE results shown, symbols correspond to the 
scaled values of the polydispersity 5 used in the simula- 
tions. Qualitatively the same behaviour is observed in 
this range. Quantitatively, the changes in volume frac- 
tion with the polydispersity 5 are smaller than in the 
soft sphere system. For larger polydispersities, the exis- 
tence of a terminal polydispersity beyond which a solid 
will be unstable manifests itself: the shadow solid curve 
never goes above 5 w 0.06 and bends down beyond this 
point. Cloud and shadow curves also begin to deviate 
from one another: the perturbative prediction that the 
curves should be identical for systems like the ones stud- 
ied here applies only to the leading 0(5 2 ) dependence of 
the volume fractions. 

Much greater differences between cloud and shadow 
curves are apparent in the number density-vs- 
polydispersity representation (Fig. 3). Here we see 
that cloud and shadow curves separate strongly as poly- 
dispersity increases. Interestingly, the corresponding 
number density can then become less than that of the 
fluid shadow. This is a signature of fractionation: the 
fluid typically contains more of the smaller particles and 
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FIG. 2. (a) Variation in the volume fraction of the cloud 
and shadow phases with their polydispersity at the freezing 
transition of soft spheres, Eq. (3), for the top- hat parent dis- 
tribution, Eq. (4). (b) MFE calculation of cloud and shadow 
curves for hard spheres with the same parent form. The sym- 
bols shown are for the scaled analogues of the values of 8 used 
in the simulations; dashed lines indicate MFE results beyond 
this range. Note that shadow phase volume fractions are plot- 
ted against their own S, not that of the parent. Dotted lines 
show exemplary tie lines connecting coexisting cloud-shadow 
pairs. 



the solid more of the larger ones as we will see below. 
As a consequence, a fluid phase can have higher density 
but lower volume fraction than a solid. 

The same reasoning explains the relative position of 
the solid shadow and solid cloud curves. A solid at the 
cloud point has the parental size distribution while a solid 
shadow phase, which coexists with a cloud point liquid, 
has a different size distribution that contains more larger 
particles. Given that both have similar volume fractions 
as we saw above, the shadow solid must therefore have a 
smaller density than the cloud point solid. In the hard 
sphere case, this effect in fact outweighs the increase in 
cloud point volume fraction with S, so that the density 
of the shadow solid decreases rather than increases as S 
increases from zero. 

Further evidence for fractionation is shown in Fig. 4, 
where we plot the cloud point pressures against the poly- 



FIG. 3. (a) Variation in the number densities of cloud and 
shadow phases with polydispersity at the freezing transition of 
soft spheres, (b) MFE calculation of cloud and shadow curves 
for hard spheres. The solid shadow curve is at lower densities 
than the solid cloud curve because the shadow solid phase 
contains more larger particles. For hard spheres the effect is 
so pronounced that the curve bends to the left. Dotted lines 
represent exemplary tie lines. 



dispersity of the parent. For a given S, fluid-solid coex- 
istence occurs in the entire range between the two pres- 
sures shown: at the lowest pressure we are at the cloud 
point coming from the fluid side, where a small amount 
of shadow solid first appears; at the highest pressure we 
are at the solid cloud point, where there is only an in- 
finitesimal amount of (shadow) liquid left. 

To understand the extent of fractionation in more de- 
tail, we consider in Figs. 5 and 6 the density distributions 
in the shadow phases. Both figures show excellent qual- 
itative agreement between the results of the soft sphere 
simulations and the hard sphere MFE calculations. In 
Fig. 5 we plot the density distributions in the shadow 
solids, i.e. the solid phases that first form from the fluid as 
we increase the parent density. We observe the trend an- 
ticipated above: the distributions are generally increas- 
ing functions of a, containing more of the larger parti- 
cles than the coexisting fluids (which have the parental 
top-hat size distribution, given that we are at the fluid 
cloud point). A further aspect becomes apparent for 
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FIG. 4. Polydispersity versus pressure at cloud points for (a) FIG. 5. Density distribution in solid shadow phases for (a) 
simulations of soft spheres, (b) MFE calculations for hard simulations of soft spheres, (b) MFE calculations for hard 
spheres. spheres. The values of the parent S used are scaled analogues 

of each other. 



larger parental S, where the shadow solid size distribu- 
tions become narrower than the parent, being suppressed 
for both the smallest and the largest particle sizes. This 
reflects the fact that in a solid it is not possible to ac- 
commodate an arbitrarily large spread of particle sizes 
without destroying the crystalline order, which is at the 
origin of the existence of a terminal polydispersity for the 
solid. 

The shadow fluid density distributions in Fig. 6 show, 
as might have been expected, that the opposite trends 
prevail in the fluid phases, i.e. the fluids that first appear 
when coming from high densities. These size distribu- 
tions are generally biased towards smaller particles; for 
more strongly polydisperse parent distributions, they be- 
come enhanced towards the smallest and largest particle 
sizes, with a minimum in the middle. Both of these can 
be understood by extrapolation from within the coex- 
istence region, where the lever rule, Eq. (2), constrains 
the density distributions in the solid and liquid phases to 
combine to the parental one. Given that we found that 
coexisting solids contain fewer small particles, and fewer 
of both extreme sizes for the larger values of 5, liquids in 
coexistence should then have more smaller particles, and 
eventually more of both the smallest and largest, exactly 



as we observe. 

The narrowing of the size distribution in the shadow 
solids can be quantified more precisely by plotting the 
polydispersity 5 of the shadow solids against the polydis- 
persity of the parent, which we here call 5^ for clarity. 
Fig. 7 shows that indeed the shadow solid has a lower 
polydispersity than the parent, an effect that gets more 
pronounced as S^ ' increases. In the MFE calculations, 
where we can reach higher parent polydispersities, the 
solid shadow polydispersity in fact starts to decrease, a 
point already observed in Fig. 2. Conversely, the accumu- 
lation of the smallest and largest particles in the shadow 
liquid causes this phase to have a larger polydispersity 
than the parent. 



VI. CONCLUSIONS 

In summary we have emphasised the necessity of fully 
catering for fractionation effects when performing the- 
oretical and computational studies of phase equilibria 
in polydisperse systems. Failure to do so can - we be- 
lieve - lead to qualitatively incorrect determinations of 
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FIG. 6. Density distribution in fluid shadow phases for (a) 
simulations of soft spheres, (b) MFE calculations for hard 
spheres. 
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FIG. 7. Polydispersities S of shadow phases against the 
polydispersity of the parent, <5^ ', for (a) simulations of soft 
spheres, (b) MFE calculations for hard spheres. 



phase behaviour. In the context of simulation studies 
we have set out in detail the reasons why unconstrained 
density ensembles such as the (/z(cr),V,T) ensemble or 
the SGCE are superior to standard (N, V, E) , (N,V,T) 
or (TV, p, T) ensembles with regard to their treatment of 
fractionation. The benefits of unconstrained ensembles 
were illustrated via an accurate determination of the 
polydispersity dependence of the freezing properties of 
size-disperse soft spheres. Our results, which are in good 
qualitative accord with moment free energy method cal- 
culations for hard spheres, show considerable fractiona- 
tion of the incipient shadow phase that coexists with the 
parent when the coexistence region is entered from the 
fluid or solid side. Moreover, we find no evidence of the 
strong narrowing of the coexistence region with increas- 
ing polydispersity that is predicted by some theories 15 ' 18 
and simulations 15 ' 26,31 ' 32 that do not account fully for 
fractionation. Although our study considered only one 
particular form of the parent distribution (a top hat), 
previous MFE studies of other parent forms (triangu- 
lar and Schulz distributions) find qualitatively similar 
behaviour. 38 We thus believe our findings for the nature 
of the freezing transition to be general. 
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